An algorithm for detecting oscillatory behavior in discretized data: the damped-oscillator oscillator detector

We present a simple algorithm for detecting oscillatory behavior in discrete data. The data is used as an input driving force acting on a set of simulated damped oscillators. By monitoring the energy of the simulated oscillators, we can detect oscillatory behavior in data. In application to in vivo deep brain basal ganglia recordings, we found sharp peaks in the spectrum at 20 and 70 Hz. The algorithm is also compared to the conventional fast Fourier transform and circular statistics techniques using computer generated model data, and is found to be comparable to or better than fast Fourier transform in test cases. Circular statistics performed poorly in our tests. Introduction Oscillatory behavior underlies many physical and biological phenomena, including the electrical behavior of neurons and neural systems. Data from neural recordings are often saved in the form of a list of discrete times when a given neuron fires an action potential. Such discretized data often appear randomly distributed in time, i.e., the spiking behavior is often well-described as a Poisson process for which the distribution of interspike intervals then is exponential. However, sometimes these data contain hidden oscillations, not obvious to casual inspection. These oscillations can be revealed by standard procedures, such as the Fast Fourier Transform (FFT) and circular statistics (CS), or by more advanced procedures including the multi-taper method [1-3]. These methods all have their strengths and weaknesses. Here we propose a simple alternative. In our approach, the experimental data are employed as an input driving force acting on a set of lightly damped computer-simulated harmonic oscillators. These oscillators we refer to as the “mathematical” oscillators. If there is a hidden oscillator in the data (a “data oscillator”), then it will drive a mathematical oscillator in a resonant way if their frequencies are the same. By monitoring the time rate of change of the energy of each mathematical oscillator, we can then detect when a data oscillator appears and disappears. The appearance of a data oscillator will cause a change in the energy of the mathematical oscillator of the same frequency, and the disappearance of that data oscillator will cause the energy of the corresponding mathematical oscillator to stop changing. Furthermore, if we lightly damp the mathematical oscillator, then we can inhibit artifactual responses of the mathematical oscillator. We refer to our method as the damped-oscillator oscillator detector (DOOD). When tested on a time series recorded from the basal ganglia of a non-human primate, the DOOD algorithm picked out oscillators at 20 and 70 Hz, while neither FFT nor CS was able to do so. We also tested the DOOD algorithm on computer generated data, to simulate a 70 Hz periodic process in the presence of Poisson noise. We found the DOOD algorithm to be superior to FFT, while FFT was superior to CS. Methods The data h(t) consist of either a continuous stream of values at every time instant t with time increment dt, or else it may consist of simply a list of times when discrete events occur. We consider the latter case for illustration. The experimental data is derived from electric potential recordings acquired from electrodes inserted into the basal ganglia of non-human primates. Whenever a unit neuronal action potential is identified, the time of discharge is automatically recorded with a precision of 0.04 msec, along with a code identifying each individual neuronal unit. The amplitude of the action potential is not recorded, nor is the duration. Thus we represent neuronal discharges by taking h(t) = h0 = 1 for a duration of 1 msec for each neuronal discharge, and h(t) = 0 otherwise. Turning our attention to the mathematical oscillators, the n th mathematical oscillator has frequency f(n) = n (in Hz), with n = 1 to N and mass taken to be unity in arbitrary units. Neglecting for now the index for n, the equation of motion for the n th mathematical oscillator is then ) ( 2 ) ( ) ( ) ( 2 0 t x g t x t h t x & & & − − = ω Eq (1) Here 0 ω is the natural frequency of oscillation for this oscillator, in the absence of friction, while g is the friction constant. We break the data down into time windows small enough such that h(t) is a constant within each time window. Consider the time window [ta, tb). The solution for Eq (1), within this time window, is then )] sin( ) cos( )[ exp( ) ( 2 0 t b t a gt h t x ω ω ω + − + = Eq (2)